clear; clc;
tic;
% 参数
N = 50;           % 中心点数
n = 100;          % 评估点数
xk = linspace(0, 1, N)';
xe = linspace(0, 1, n)';
F = @(x) sin(5*pi*x) .* exp(-2*x);
y = F(xe);

M = 20;           % 最大基函数个数
loo_history = [];

%% LOOCV优化epsilon
epsilon_list = linspace(0.01, 20, 100); % 候选epsilon范围
cv_errors = zeros(size(epsilon_list));
% 预先计算距离矩阵（n x N）
DM = pdist2(xe, xk); % 或者使用abs(xe - xk')，但需要确认维度
for i = 1:length(epsilon_list)
    epsilon = epsilon_list(i);
    Phi = exp(-epsilon * DM.^2); % 生成基函数矩阵
    
    % 计算系数
    w = pinv(Phi) * y;
    
    % 计算帽子矩阵的对角元素
    H = Phi * pinv(Phi);
    H_ii = diag(H);
    
    % LOOCV残差
    residual_loo = (y - Phi * w) ./ (1 - H_ii);
    cv_errors(i) = mean(residual_loo.^2);
end
% 选择最优epsilon
[~, idx] = min(cv_errors);
epsilon_opt = epsilon_list(idx);
fprintf('LOOCV优化结果：最优epsilon=%.4f\n', epsilon_opt);
% 使用最优epsilon重新生成Phi矩阵
epsilon = epsilon_opt;
Phi = exp(-epsilon * DM.^2);  

%% 
% 基函数矩阵
Phi = exp(-epsilon * (pdist2(xe, xk)).^2);

%% 贪心算法初始化
I = [];     % 已选中心索引
A = [];     % 设计矩阵
C = [];     % 中心点坐标
r = y;      % 初始残差

%% 带LOOCV的贪心迭代
for k = 1:M
    available = setdiff(1:N, I);
    if isempty(available)
        break; 
    end
    
    % 选择投影最大的基函数
    inner_products = Phi(:, available)' * r;
    [~, idx] = max(abs(inner_products));
    best = available(idx);
    
    % 更新设计矩阵和中心集合
    I = [I, best];
    A = [A, Phi(:, best)];
    C = [C; xk(best)];
    
    % QR分解求解权重
    [Q, R] = qr(A, 0);
    w = R \ (Q' * y);
    r = y - A * w;  % 更新残差
    
    % 计算当前模型的LOOCV误差
    H_ii = sum(Q.^2, 2);
    loo_e = r ./ (1 - H_ii);
    loo_mse = mean(loo_e.^2);
    loo_history = [loo_history, loo_mse];
    
    % 终止条件：LOOCV误差上升
    if length(loo_history) >= 2 && loo_history(end) > loo_history(end-1)
        fprintf('LOOCV误差上升，终止于迭代 %d\n', k);
        % 回退到前一步
        I = I(1:end-1);
        A = A(:, 1:end-1);
        C = C(1:end-1);
        break;
    end
end

%% 最终模型（回退后）
[Q, R] = qr(A, 0);
w = R \ (Q' * y);
fe = A * w;

%% 可视化
yy_true = F(xe);
yy_pred = fe;



%  LOOCV误差曲线
figure;
semilogy(epsilon_list, cv_errors, 'b-', 'LineWidth', 1.5);
xlabel('形状参数γ'); ylabel('LOOCV误差');
title('LOOCV误差随γ变化曲线');
grid off;
figure;

plot(xe, yy_true, 'k-', 'LineWidth', 1);
hold on;
plot(xe, yy_pred, 'r--', 'LineWidth', 1.5);
xlabel('x'); ylabel('y');
legend('真实函数', '稀疏RBF近似');
title('使用LOOCV的贪心算法稀疏化RBF近似');
grid on;

% 误差计算
error = max(abs(yy_true - yy_pred));
rmse=sqrt(mean((yy_true - yy_pred).^2));         %均方根误差
fprintf('最大绝对误差: %.2e\n', error);
fprintf('RMSE: %.2e\n', rmse);
disp(['使用的基函数数量: ', num2str(length(C))]);
toc;